Bridging Time Scales in Cellular Decision 
Making with a Stochastic Bistable Switch 

Steffen Waldherr, Jingbo Wu, and Frank Allgower 

Institute for Systems Theory and Automatic Control, 

Universitat Stuttgart, PfafFenwaldring 9, 

Stuttgart, Germany 

f— s Cellular transformations which involve a significant phenotypical change of the cell's 

,— I state use bistable biochemical switches as underlying decision systems. 

^^ In this work, we aim at linking cellular decisions taking place on a time scale of years 

CN to decades with the biochemical dynamics in signal transduction and gene regulation, 

^-^ occuring on a time scale of minutes to hours. We show that a stochastic bistable switch 

^ C^ forms a viable biochemical mechanism to implement decision processes on long time 

^> scales. As a case study, the mechanism is applied to model the initiation of follicle growth 

_ in mammalian ovaries, where the physiological time scale of follicle pool depletion is on 

^•^ the order of the organism's lifespan. We construct a simple mathematical model for this 

process based on experimental evidence for the involved genetic mechanisms. 

I""! Despite the underlying stochasticity, the proposed mechanism turns out to yield re- 

PP liable behavior in large populations of cells subject to the considered decision process. 

\^ Our model explains how the physiological time constant may emerge from the intrinsic 

Q stochasticity of the underlying gene regulatory network. Apart from ovarian follicles, the 

proposed mechanism may also be of relevance for other physiological systems where cells 

take binary decisions over a long time scale. 






^ 



Background 



> 

j/^ The dynamics of biological systems span a wide range of temporal and spatial scales. The interactions 

\^ among dynamical properties on different scales govern the overall behavior of the biological system, 

^^ and thus form an important area of computational research in biology |18j . A particularly interesting 

^~^ question in this context is how the behavior on a slow time scale emerges mechanistically from the 

l/~j dynamics on fast time scales. For example, how do cell population dynamics in tissues, which may 

^D evolve on a time scale of months, years or even decades, originate from the dynamics of the underlying 

^—•^ gene regulatory networks, with a time scale of just minutes to hours? 

. . In this work, we aim at bridging the time scale from gene regulation to cellular transformation 

>>" processes on the tissue or cell population level. We specifically consider cellular transformation 

processes based on a bistable biochemical switch. Such switches have two distinct stable stationary 
states, and the cell initiates a transformation when the switch changes from one stable state to 
Cd the other one. Bistable switches have previously been used to model a large number of cellular 

transformation events, such as progression through cell cycle arrest in the maturation of Xenopus 
oocytes [TTl [TU] or initiation of programmed cell death [HI and cellular differentiation [3] in higher 
organisms. Most models for these systems are constructed as deterministic models, and thus an 
external stimulus is required to induce changes in the switch's state. 

Here, we consider transformation processes which apparently do not require any external stimulus 
to be initiated, and which still follow reliable temporal characteristics. Reliable thereby means that 



in a large population of cells, the number of cells that have already initiated the transformation 
can be described deterministically with high accuracy. We propose a generic transformation process, 
where a phenotypical change in the state of a cell is initiated as soon as a bistable biochemical switch 
changes its internal state. In previous studies, random switching caused by internal fluctuations is 
usually attributed to pathological events |15j . In the mechanism proposed here, random switching 
has a regular physiological function. 

A striking example for the kind of transformation processes we aim to describe is involved in mam- 
malian oocyte maturation. In mammalian females, all or almost all of the oocytes that will ovulate 
through the organism's life-span are already present at birth or shortly thereafter as a population 
of so-called primordial follicles. Throughout the organism's reproductive life, follicles undergo the 
primordial to primary transition, which marks the start of a development process that will eventually 
lead to either ovulation or removal of the oocyte through atresia [12 [5^. In this way, there is a steady 
supply of mature follicles for ovulation, while the pool of primordial follicles is gradually depleted. 
The mechanisms through which the follicle transition is initiated are largely unknown, although a 
number of ovarian factors that may be relevant have been identified experimentally [231 [H [3T] . Im- 
portantly, the transition seems to be regulated locally in the ovary, and not through the endocrine 
system ^. An astonishing observation in this process is that in one follicle, the transition may occur 
already a few months after generation of the primordial follicle pool, while another follicle may stay 
several decades (for organisms with a sufficiently long lifespan) in the resting stage before growth is 
initiated. From the medical side, a misregulation of this process is implicated in premature ovarian 
failure due to follicle depletion, which is a major reason for infertility in human females. 

By way of a case study, we apply the proposed transformation mechanism to the problem of 
growth initiation in ovarian follicles. Including also cell-cell interactions supported by experimental 
evidence, we obtain a physiologically plausible model for this process, showing very good agreement 
with human clinical data on a scale of several decades. 

Models and Methods 
Deterministic model of a bistable switch 

The model of a bistable switch that we use is based on a positive feedback loop between two com- 
ponents. Consider a biochemical reaction network involving the two molecular species X and Y. 
Mathematically, the temporal evolution of the amounts of the two species is described with the 
ordinary differential equation 



i = vi+ V2{y) - V3{x) ^ ki + j^h , h " "1^ 
y = v^ix) - v^{y) = ~U2y, 



(1) 



where x and y denote the amounts of X and Y, respectively. The vector (a;, y) will be referred to 
as the microstate of the biochemical reaction system. Ultrasensitivity, which is required to achieve 
bistability |11| . is generated by the Hill- type production rates V2 and W4. 

In the sequel, we will assume that the molecular species X and Y represent gene transcripts, and 
the amounts x and y indicate the respective transcript copy number. The nominal parameter values 
that we use are given in Table [l] For simplicity, we assume that the parameters are symmetric, 
i.e. Vi = V2, Ml — M2 and ui — U2. The parameter values are within the physiological range 
for typical gene transcription processes. In particular, the degradation rate of 0.01 ^^ corresponds 
to a gene transcript half-life time of about 70 minutes. Typical transcript half-life times are in a 
range from several minutes to several hours. The minimal transcription rate of X is given by ki and 
corresponds to 3.3 transcripts that are produced per hour. The transcription rate upon maximal 



Table 1: Nominal parameter values for the bistable switch model (IT]). Transcript copy numbers are 
considered to be dimensionless. 



Parameter Value 



mm 



1^1 0:05^ 

Ml 2 25 

0-01 ± 



Parameter Value 



^^2 0:55^ 

h 3 



mm 



activation is given by Vi_2 and corresponds to 33 transcripts produced per hour. This is an arguably 
low transcription rate, but it was mainly chosen to allow for an efficient computational treatment in 
the stochastic case. The basic principle also works with higher maximal transcription rates. 

For two-dimensional systems, it is convenient to check bistability by considering nullclines in the 
state space [7]. With this graphical representation, it is also easy to evaluate how good the two stable 
states are actually separated [3] . The nullclines for the model given in M , with nominal parameter 
values, are depicted in Figure [TK. From the figure, it is clear that there are three equilibrium points, 
labelled I, II and III. A stability analysis of the equilibrium points shows that the deterministic 
system described by (IT]) is bistable, and the corresponding reaction network implements a bistable 
switch. We construct a macrostate for this system by defining the three sets ^off , ^on , ^trans C M^ , 
corresponding to the switch being ojf, on, or in transition, respectively, ^off contains the equilibrium 
point I, Qtrans contaius II, and flon contains III. Moreover, to have a well-defined macrostate for 
each microstate that the system can reach, we require that ^off U flon U fltrans — H^+i where M+ 
stands for the non-negative real numbers. For our model, we define 

noff^{{x,y)eRl \x + y<k} 
non^{{x,y) eM.l\x + y>h} 

and, consequently, 

^tra,is ^ {{x,y) e r1 \ h < X + y < 1.2}, 

with suitable parameters li and l2- With model parameters as given in Table IT] a suitable choice 
which we will use in this work is Zi = 25 and ^2 = 55. 

Stochastic model of a bistable switch 

The deterministic model of the bistable switch discussed in the previous section is suitable to describe 
the existence of two distinct macrostates, corresponding to stable equilibrium points in the model. 
However, to capture transitions between these macrostates which are caused by intrinsic fluctuations, 
a stochastic model has to be considered. The stochastic description of the bistable switch makes use 
of the Markov property of biochemical reaction networks. In a stochastic setting, the amounts of 
molecular species may only take discrete values from the set 17 = {{x,y)'^ \ x € Nq,?/ S Nq}. The 
stochastic state of the switch at time t is given by the discrete probability distribution P{x,y,t), 
which for each microstate {x,y) € 17 gives the probability that the switch is in this microstate at time 
t. To describe the temporal evolution of the probability distribution, we use the chemical master 
equation (CME) [28]. The reaction network for the bistable switch is not described with elementary 
reactions only, and thus it is not possible to construct the CME according to its rigorous derivation 
|13j . However, a theoretical investigation by Rao and Arkin |24j has shown that as an approximation, 
the propensity functions for state transitions can be taken from the according reaction rate laws. 





•c + y = h 



Figure 1: Characterisation of the phase space in the bistable switch model (IT]). A: Nullclines for the 
deterministic model of the bistable switch. / and /// are stable equilibrium points, // is an 
unstable one. B: Schematic illustration of the configuration space for the Markov process 
Q describing the cell transformation process. 

For the bistable switch described above, we thus use the CME 
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for {x, y) £ il. 

In the stochastic description, we can compute the probabilities that the switch is in any of its 
three macrostates directly from a solution of the CME. Define Poff{t), ptrans{t) and Pon{t) as the 
probabilities that the switch is off, in transition and on, respectively. Given a solution of the CME, 
these can be computed by summing up the probabilities that the system is in the corresponding 
microstates, i.e. Pon{t) = ^{x i/)en ^(^j2/j*)j ^^d equivalently for the other macrostates. 



A transformation process modelled with a stochastic switch 

Cellular transformation processes are often based on a bistable biochemical or genetic switch. In 
the initial state of the cell, the switch would be in the off state. Switching to the on state implies 
a significant change in the amount of an involved signaling molecule, e.g. a transcription factor. If 
the on state is maintained for some time, this change would result in a larger phenotypical change 
of the cell, e.g. through significant changes in gene expression. The mechanisms that induce this 
change are not part of the stochastic switch, but from a signaling perspective downstream of it. 

Most transformation processes rely on specific external stimuli, and the cell will initiate the trans- 
formation upon encountering the required stimulus. There are however examples where such a 
stimulus is not strictly required, and this is the case that we are dealing with in this paper. More- 
over, we will focus on the behavior of cell populations, studying the problem how the temporal 
dynamics of the transformation process evolve in a pool of many cells. 

The basic mechanism that actually triggers the bistable switch in our model without an external 
stimulus are the intrinsic fiuctuations of concentrations in any biochemical reaction network, that 
are due to the stochastic nature of chemical reactions. As a rare event, these fiuctuations may 
become so large that the microstate of the system crosses the separatrix between the domains of 
attraction in the deterministic system. As a consequence, the microstate around the other stable 



equilibrium point will become strongly attractive, and the switch will change its niacrostate to on 
with a high probability. In this paper, we assume that the transformation is irreversible, which fits 
well to the process of follicle growth initiation. Also other processes such as programmed cell death 
are irreversible. 

The described transformation process is easily modelled as a continuous-time Markov process. If 
the switch is in one of the macrostates ojf or in transition, then we directly use the microstates and 
transition probabilities of the underlying biochemical reaction network to model the transformation 
process. To account for the irreversibility of the transformation, the on state of the switch corre- 
sponds to only one state of the Markov process, which is an absorbing state. The transitions of other 
microstates to the absorbing state are governed by the propensity functions for the corresponding 
transitions in the underlying biochemical network. The resulting state space for the Markov process 
model of the transformation process is shown in Figure [lj3. 

In our model of the stochastic switch, the macrostates off and in transition are defined by compact 
regions in state space. As a consequence, the Markov model of the considered transformation process 
has a finite state space, and can therefore be treated computationally with standard approaches. Let 
P(x,y,t) G M" denote the complete probability state vector of the system with n discrete states 

(x,y), 

P(x,y,<) = (P(0,0,t), F(l,0,i), P(0,l,t), ..., -P( — ,0,t), Pixc,yc,t)f, (3) 

Ui 

where (xclJc) stands for the microstates which correspond to the on state of the switch (see Fig- 
ure [l|3). The on state is made irreversible by removing the associated transitions from the Markov 
process. The master equation can be written as the linear ordinary differential equation 

P(x,y,i) = AP(x,y,t), (4) 

where A G M"^" is the state transition matrix. The matrix A can be computed directly from the 
values of the reaction propensity functions in each microstate |19| . 

The differential equation Q can be solved using standard tools for numerical integration. For the 
results described in this paper, we used the odelSs function in Matlab (The MathWorks, Natick, 
MA) to obtain a numerical solution of Q. 

Results and Discussion 

A hypothetical mechanism for oocyte maturation 

In this section, we suggest a biochemical mechanism that offers a molecular explanation for the large 
depletion times of several decades in the human oocyte pool. The model is based on experimental 
evidence obtained in a very informative series of studies by Skinner and colleagues (see |26| for a 
review), where the influence of several ovarian factors on the primordial to primary transition as well 
as some interactions between them have been studied. Because a positive feedback loop is necessary 
for a bistable switch [T6^, we have specifically searched for such an interconnection. 

Indeed, experimental evidence suggests a positive feedback circuit involving two ovarian factors 
that are relevant in the primordial to primary transition: the factor KIT ligand (KITL) is produced 
by granulosa cells surrounding the oocyte and stimulates both the oocyte and surrounding theca cells 
to promote follicle development. Moreover, KITL stimulates the production of both keratinocyte 
growth factor (KGF) and hepatocyte growth factor (HGF) in the surrounding theca cells. KGF 
and HGF themself stimulate the production of KITL in the granulosa cells, thus providing a posi- 
tive feedback loop 22 . Moreover, the oocyte of primordial and developing follicles produces basic 
fibroblast growth factor (bFGF), which acts on surrounding granulosa cells and has been shown to 
increase the expression of KITL [5T] . 




Figure 2; Hypothetical biochemical network for the primordial to primary transition in ovarian 
follicles 

These pieces of experimental evidence thus support the hypothetical mechanism that is shown 
in Figure ^ Our simplistic mathematical model presented in ([l]) can be used to describe this 
mechanism, where the variable x represents granulosa-derived KITL activity and y represents theca- 
derived KGF and HGF activity. The reaction vi describes the influence on KITL expression of 
oocyte-derived bFGF, which is here assumed to be constant. The reactions V2 and U4 arise from the 
positive feedback interconnection, whereas v^ and v^ describe a constitutive degradation of KITL, 
KGF and HGF. 

The stochastic switch generates reliable long-term behavior 

The differential equation Q that governs the initiation probabilities of the irreversible transformation 
process is a linear ordinary differential equation, so in principle it can be solved analytically. Due to 
the size of the system, this is however not feasible. Yet, we can characterize the probability that a 
given cell has initiated the transformation process by the explicit formula 

n 

Pon{t) = 1 - cie-^i* + Y, c^{t)e-^^\ (5) 

4 = 2 

where Ci > 0, < Ai < Re(Ai) for i — 2, . . . ,n, and the Ci(t) are polynomials in t. The mathematical 
derivation of ^ is provided in the appendix. 

From considering the general form of the analytical solution given in ([5]) , we obtain two important 
conclusions about the stochastic transformation process. First, we observe that the probability for 
a given cell to initiate the transformation tends to 1 as time increases. Second, because Ai is the 
dominant decay rate, for larger times t ^ the probability of not having initiated the transformation 
can be approximated by 1 — Pon{t) ~ Cie""^^*, a simple exponential decay. For the biochemical 
parameter values given in Table [11 the numerical solution for Pon{t) is shown in Figure [sK. For these 
parameter values, which are in the physiological range for the considered biological processes, we 
indeed get to a time scale of years to decades in the probability of the transformation event, with a 
half-life time of about 5.9 years. 

Let us now move to the population level, and consider a pool of cells, each of them being subject to 
the considered transformation process with a bistable switch. In the first step, we make the simplistic 
assumption that no interactions among the cells are taking place, so individual transformations are 
probabilistically independent events. The number of remaining cells Nr(t) can easily be characterized 
by a binomial distribution as 

P{NAt) ^N)= (J^^ (1 - Pon{t)fpon{t)'''>-''. (6) 
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Figure 3: Dynamical characteristics of the stochastic bistable switch on the single cell level and 
the population level. A: Probability of transformation event Pon{t)- B: Population size 
probability distribution over time. C: Probability density function of the depletion time 



The properties of the binomial distribution give the expected number of ceUs remaining in the pool 
as 

E[Nrit)]=Noil^Pon{t)), (7) 

where Nq is the initial number of cells in the pool. 

The probability distribution P{Nr{t) = N) for the population size in the transformation process 
considered in this paper is shown in Figure l3j3 as a function of both cell number N and time t. The 
number of initial cells Nq = 10^ was chosen from the reported range of ovarian follicles, 7 • 10^ to 
2 • 10^, in human females at birth. For each point in time, the distribution has a very sharp peak, 
which indicates that the average value E[Nr{t)] is a reliable prediction for the number of cells that 
have already undergone the transformation at a given time. 

A relevant characteristic of the considered process is the time at which the initial cell population is 
depleted, i.e. when nearly all cells have undergone the transformation. To make this notion precise, 
we introduce the depletion number Nd- The depletion time Td is defined as the smallest time t such 
that Nr{t) < Nd, i.e. only Nd cells are remaining in the initial population. For the process of follicle 
growth initiation, we use Nd = 1,000, which has been considered to mark the onset of menopause 

i- 

The cumulative probability distribution function for the depletion time Td is computed from the 
distribution obtained in Q as 

P{Td<t)=P{Nrit)<Nd)= J2 PiNrit)^N). ^g) 

N<No 

The probability density function for the depletion time is computed by taking the derivative of 
the cumulative probability distribution function Q. The resulting probability density function 
for the depletion time in follicle growth initiation is shown in Figure |3p. From the density func- 
tion, the expected value and the standard deviation are obtained as ^'[T'd] — 59.1 years and 
a/£;[T|] - E[Td]^ = 0.27 years, respectively. 

The expected value for Td can also be computed by solving 1 — Pon{Td) — jf- Using ([s]), it can 
thus be approximated by E[Td\ ~ 3^ In -^, where Ai is the dominant decay rate of the process. 

Next, we compare the computed statistical characteristics of the follicle depletion process to med- 
ical data. Explicit follicle counts are only sparsely available. However, the available pieces of data 
indicate that fluctuations in actual follicle numbers are larger than predicted by our model |14j . 
Concerning the depletion time, a recent medical study suggest an average age of 51.1 years for the 
onset of menopause, with a standard deviation of 3.8 years [5]. Our model predicts a depletion 
time of Td — 59.1 years, which is reasonably close to the experimentally observed depletion time. 
However, the standard deviation of 0.27 years in our model is significantly less than observed from 
medical data. In summary, even though our model is based on a highly stochastic process, the anal- 
ysis reveals that it leads to much more reliable temporal characteristics than observed in the real 
system. This indicates that stochastic effects alone may not be sufficient to explain the heterogeneity 
observed in the follicle depletion process. 

An alternative explanation would be by heterogeneous parameter values among individual or- 
ganisms. This explanation is also supported by statistical analyses of medical data [5 , where it is 
suggested that the onset of menopause is largely based on genetic factors, which would be related 
to parameter values in our model. To investigate this possibility, we have computed the expected 
depletion times for different parameter values. The computation was based on the eigenvalues _of 
the transition matrix A and the approximation £'[7^] « 3^ In j^. The results are given in Table 
From these results, we note that even small parameter variations in the model of the bistable switc 
lead to very large variations in the expected depletion time. This is not realistic for a biological sys- 
tem, and in the following section we explore mechanisms to increase the robustness of the depletion 
time with respect to parameter variations. 



Table 2: Expected depletion times (years) in the single cell model Q. Expected depletion times 
(years) in the single cell model m for multiplicative variations in single parameters. For 
simplicity, we always assume Vi = V2, ui = U2, and Mi = M2- 



Param. 


Factor 


0.8 


0.9 
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1.1 


1.3 


fci 
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Vl,2 
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0.4 
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1010 
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Mi,2 




0.2 


2.1 


10 


418 


3410 


2.0-10^ 


h 




0.6 


6.3 


20 


163 


420 


9.6 • 103 



Increased robustness by interactions on the population level 

In the last section, we have characterized the properties of the transformation process based on a 
bistable switch, with the depletion time of a pool of cells subject to the transformation as character- 
istic output of the model. We have shown that the proposed model produces reliable depletion times, 
in the sense of a small standard deviation, for fixed values of the biochemical parameters. However, 
we have also observed that the average depletion time in the basic model is quite sensitive to vari- 
ations in the biochemical parameters. Clearly, this large sensitivity is not acceptable in a model 
that should be a meaningful representation of the primordial to primary follicle transition. In this 
section, we propose an additional mechanism that reduces the sensitivity of the average depletion 
time significantly. 

The additional mechanism is based on the experimental observation that follicles in later stages 
of development actively suppress the primordial to primary transition in resting follicles '26] . The 
inhibition of follicle growth initiation is mediated by the Anti-Miillerian hormone (AMH), which 
interferes with stimulatory signals by KITL, bFGF, and KGF [20 . Although AMH is known to 
signal via SMAD proteins |5S| , the molecular mechanisms of follicle growth inhibition by AMH seem 
to be unknown. To include the inhibitory effect into the simplistic switch model (IT]), we assume 
that the rate of KITL production in primordial follicles is reduced with an increasing number of 
growing follicles. This is achieved by changing fci in the original model given in (fTj) from a constant 
parameter to the expression 

ki{n2) = -r^ , (9) 

A„ + 712 

where ki^max is the maximal production rate of KITL, n2 is the number of growing follicles, and Kn 
is an additional parameter. While follicle development is a complex process with many intermediate 
stages |30j . in this analysis we use a simple two-state population model, where ni denotes the number 
of primordial follicles, and n2 the number of growing follicles. The assumptions of the model are 
that primordial follicles initiate growth with a rate as determined by Ai in (O. Growing follicles are 
assumed to stay in this stage for a constant amount of time r, after which they leave the pool either 
through ovulation or atresia. From these specifications, one can derive a model given by the system 
of delay-differential equations 

ni(t) = -Ai(n2(t))ni(t) 

n2[i) = Ai(7i2(t))ni(t) - Ai(n2(i - T))ni(t - r), 

where Ai(n2) is the decay rate computed from the transition matrix A in Q. 

Using the parameters in TablelSJ the population model given by ( 10 ) now predicts a depletion time 



of Td = 50.0 years, which is almost equal to the depletion time suggested by the medical study [S]. 



The development of the ovarian follicle pool over time, as predicted by the model in (10 1, is shown 



in Figure [4] The prediction is compared to clinical data of follicle numbers at different ages taken 



Table 3: Nominal parameter values for the population model (10) 
Parameter Value Parameter Value 
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Figure 4: Evolution of follicle number. Model predictions from (10) (line) vs. clinical data from [8] 
(crosses). 



from [8]. Although the parameters have only been adjusted to fit the depletion time, the predicted 
time course is reasonable close to the clinical data. In particular, the proposed population model 
(fTol) intrinsically captures the previously observed increase in the follicle depletion rate at an age of 
approximatively 38 years [S]. 

In order to investigate the sensitivity of the extended model to variations in the biochemical 
parameters, we have computed again the expected depletion times for different parameter values. 
The results are given in Table |4J The variation in the depletion time is significantly reduced, 
compared to the model ^, where follicle interactions are neglected. It should also be pointed out 
that the depletion time is quite insensitive towards variations in the two parameters K^ and t 
which were newly introduced in the population model. This result illustrates that the robustness of 
the depletion time with respect to parameter variations may be substantially increased by adding 
interactions among individual follicles to the proposed model of the transformation process. 



Conclusions 

In this paper, we deal with a fundamental question in the development of multicellular organisms: 
How can biochemical reactions and genetic mechanism acting on the scale of minutes in individual 
cells generate dynamics with characteristic times of years to decades on the tissue level? As a possible 
mechanism to achieve this, we propose a generic transformation process based on a bistable stochas- 
tic switch. From the underlying genetic interactions and biochemical reactions, the process can be 
modelled as a continuous-time Markov process. We show that the proposed stochastic mechanism 
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Table 4: Expected depletion times (years) in the population model ( 10 ) for multiplicative variations 
in single parameters. 



Param. F^^tor 
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generates reliable long-term behavior on the population level, with cells undergoing the transforma- 
tion with an exponentially decaying rate. Thereby, the decay rate is equal to the dominant eigenvalue 
of the transition matrix describing the underlying biochemical network. 

We pose the hypothesis that a biological instance of this mechanism is present in the development 
of ovarian follicles. To describe this process, we constructed a simple model of a bistable switch 
in the primordial to primary transition for ovarian follicles. The model is based on experimentally 
determined factors and their interactions in the different cell types constituting the ovarian follicles. 
Although it is not assured that a bistable switch in ovarian follicles will indeed be based on the factors 
that we have used here, the basic mechanism would work equivalently well with other factors. 

Despite its simplicity, our model explains well how the long-term characteristics of follicle devel- 
opment may reliably be generated by biochemical reactions occurring on much shorter time scales. 
Keeping the model simple serves two purposes: first, it shows that the dynamics of follicle growth 
initiation can be generated by a quite simple mechanism. Clearly, additional pathways and regula- 
tory feedback interactions that we have not included in this model can be expected to be present in 
the system. These may serve to increase robustness of the network, or to provide additional inputs 
to control the transition rate, e.g. for the endocrine system. Second, the simplicity of the model 
allows us to solve the chemical master equation for the network numerically, and thus to obtain a 
quantitative description of the model behavior. 

As a possible shortcoming of the basic model on the single cell level, we observe an unrealistic 
large sensitivity of the follicle depletion time with respect to parameter variations. By adding the 
experimentally supported inhibition of follicle growth initiation by later-stage growing follicles, the 
sensitivity of the depletion time could be reduced significantly. Apart from the inhibition included 
in the model, other interactions among individual follicles seem to play a role in the primordial to 
primary transition |25j . We envision that the inclusion of more regulatory interactions may further 
decrease the sensitivity of the depletion time with respect to parameter variations to a physiologically 
plausible level. 



Appendix: Computation of the transition probability 

In this section, we prove that the probability that a given cell has undergone the considered trans- 
formation process is given by Ponit) as in ^. The proof is based on considering the solution of the 
underlying CME Q. 

Since the last microstate is an absorbing state of the Markov process, Q can be written as 



^abs 



(11) 
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where Ared G m("^i)^("^i) describes the interactions among the non-absorbing states, and aabs G 
]]jix(n-i) cj^ggcribes the transition propensities to the absorbing state. 

Let us first derive some essential properties of the matrix Ared- Since A is a stochastic matrix, we 
have 

n-l 

J2 \^J^\<~A^^, (12) 

for J = 1, . . . , n — 1, i.e. Ared is diagonally dominant. Thus, Gersgorin's theorem [T7] asserts that all 
eigenvalues of Ared have a non-positive real part. Even more, since aabs is non-zero, ( |12| holds with 
a strict inequality for at least one i. Thus, by Theorem 10.7.2 in [17], all eigenvalues of Ared have 
negative real part. By the properties of the considered biochemical network, Ared is irreducible, and 
its off-diagonal elements are non-negative. From Corollary 4.3.2 in [37], it follows that Ared has an 
eigenvalue Ai G M with algebraic multiplicity 1 and a strictly positive corresponding eigenvector vi 
such that Re A < Ai for all A 7^ Ai in the spectrum of Ared- 

Denoting Pred — {Pi, - - - , Pn-i)"^ , we have Pred = AredPred- From the previously derived prop- 
erties of the matrix Ared, the general solution of this differential equation is given by 

s 

Predit) = awie^i* + Y^ vfc,{t)e^^\ (13) 

where Ci{t) are polynomials in t and a is a constant coefhcient, depending on the initial condition 
Pred(O). The condition Pred{t) > for all t implies that a > 0. The transition probability Pon{t) is 
computed as 

Ponit) = 1 - 1^ Predit) 

= l-ae^^' + Y,c^{t)e^•', ^^^^ 

i=2 

where a = al"""!;! > and Ci{t) = Ci{t)l^Vi, 1 = (1, . . . , 1)""". 
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